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Abstract. - Dissipative particle dynamics (DPD) and multi-particle collision (MPC) dynam- 
ics are powerful tools to study mesoscale hydrodynamic phenomena accompanied by thermal 
fluctuations. To understand the advantages of these types of mesoscale simulation techniques in 
more detail, we propose new two methods, which are intermediate between DPD and MPC — 
DPD with a multibody thermostat (DPD-MT), and MPC-Langevin dynamics (MPC-LD). The 
key features are applying a Langevin thermostat to the relative velocities of pairs of particles or 
multi-particle collisions, and whether or not to employ collision cells. The viscosity of MPC-LD 
is derived analytically, in very good agreement with the results of numerical simulations. 



Soft matter systems such as polymer solutions, colloidal suspensions, vesicles, cells, and 
microemulsions exhibit many interesting dynamical behaviors, where hydrodynamic flow plays 
an important role, as do thermal fluctuations. Several mesoscale simulation techniques for 
the flow of complex fluids accompanied by thermal fluctuations have been developed in the 
last decades, such as direct simulation Monte Carlo (DSMC) [1,2], the Lattice Boltzmann 
method [3], dissipative particle dynamics (DPD) [4-12], and multi-particle collision (MPC) 
dynamics [13-17]. These methods have many similarities. The most important common 
feature is the local mass and momentum conservation, which is crucial to obtain hydrodynamic 
behavior in the continuum limit. Since most of these methods were developed independently, 
the relations between them are not well explored so far. In this letter, we clarify the relations 
between the particle-based off-lattice methods, particularly DPD and MPC, and use this 
insight to propose two new intermediate methods (all summarized in Fig. ^| . 

We start from the relations between the Langevin and the Andersen's thermostat (AT) [18] . 
The underdamped Langevin equation of N particles is given by 

m-^- = -ViU + f LT , (1) 

/ LT = -7V( + <7&(«), (2) 

where V, = d/dr^ m is the mass of a fluid particle, and r, and Vj are the position and velocity 
of the i-th particle, respectively. The force /lt represent the Langevin thermostat. To satisfy 
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the fluctuation-dissipation theorem, the Gaussian white noise £j(t) has to have the average 
(£i,a(t)) — and the variance (£i,a(t)£j,/3(t')) = 2k^T8ij8 a p5[t — t') 1 where a,fi € {x,y,z}and 
k^T is the thermal energy, and its amplitude er in Eq. Q is related to the friction constant 
7 by 7 = (T 2 . Each particle has one thermostat (heat bath), which is independent of all 
other particles. This thermostat does not conserve momentum, and hence the hydrodynamic 
interactions are not taken into account. We separately integrate the potential forces —ViU 
and the Langevin thermostat /lt using a multiple-time-step algorithm [19], where a shorter 
time step is employed for —VJJ. The Langevin thermostat is integrated by the leapfrog 
algorithm [20], which implies 

Vj(t n +i) = avj(t n ) +6£i, n , (3) 

1 - jM/2m yfiEt/m 

where a = — ; , b = : — ; (4) 

l + 7Ai/2m l + 7Ai/2m 



with Vi(i n+1 / 2 ) = {vi(t„+i) + Vi(t„)}/2 and £i(t„+i/ 2 ) = £i,n/V At. The modified Verlet al- 
gorithm in Ref. [10] also gives Eq. 0. This thermostat works even with large time steps At > 
2m/ '7. For the thermodynamically ideal gas (J7 = 0), Eq. © gives v,(t n ) = a '&£i,r»— 
This velocity exhibits a Maxwell-Boltzmann distribution with (vi ta (t n )vi_ a (t n )) — k^T/m for 
any At. Thus, this discretized thermostat belongs to the class of generalized ATs described 
in Ref. [12]. For jAt/2m — 1, the first term in Eq. © vanishes (a = 0), and a new velocity 
v i(t n +i) is selected from a Maxwell-Boltzmann distribution. This corresponds to the origi- 
nal Andersen's thermostat [18]. Thus, AT can be interpreted as the discrete version of the 
Langevin thermostat with jAt/2m = 1. This relation between the Langevin and Andersen's 
thermostat remains valid in DPD and MPC, as shown below. 

The DPD thermostat is a modified Langevin thermostat, which applies the relative ve- 
locities of the neighbor pairs. This implies that the thermal force in Eq. is now given 
by 



/dto = I -w(nj)(vi - Vj) ■ fy + \fw(nj)^j(t) \ : 



(5) 



where = — r_y, r,j = |r,j|, = Tij/rij, and £ji(t) = — £ij(t). This thermostat is applied 
only in the direction f y to conserve the local angular momentum. In DPD, a linear weight 
function, with w(rij) = 7(1 — rij/r cut ) for ry /r cut < 1 and w(rij) = otherwise, is typically 
employed. Furthermore, DPD is usually combined with soft repulsive potentials U(rij) — 
{A/2)(\ ~ rij/r cut ) 2 with cutoff at r cut [5], but other pairwise or multibody potentials are also 
available [8]. DPD shares many features with smooth-particle hydrodynamics (SPH) [21], a 
method to solve the Navier-Stokes equation in a Lagrangian representation, and a modified 
version [9] of DPD corresponds to SPH with thermal fluctuations. In Shardlow's splitting 
algorithm [10], each thermostat of the ij pair is separately integrated. This algorithm with 
multiple time steps gives the Andersen-thermostat version of DPD proposed by Lowe [22]. 
An energy-conservation version of DPD (micro-canonical ensemble), called DPD+e, has also 
been proposed [6,7]. DPD+c needs an additional variable in order to exchange the momenta 
of particles, since five of six degrees of freedom of a particle-pair are fixed by the conservation 
of the translational (3) and angular (2) momenta. We modify the DPD thermostat to remove 
angular-momentum conservation (DPD-a) for the discussions below. In this case, the thermal 
force reads 

/dti = -w(nj ) (vj - Vj ) + ijw(rij)tij (t) . (6) 
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Fig. 1 - (Color online) Relations between particle-based hydrodynamic methods. 



We call the versions of methods with or without angular-momentum conservation '+a' or 
'-a', respectively, /dti still keeps the translational momentum conservation. Note that the 
numbers of the thermostats are NN n h/2 and 3NNnh/2 for /dto and /dti, respectively, where 
iVnb is the mean number of the neighbors with m < r CU f These can be much more than the 
number of degrees of freedom 3 A. 

MPC is a modification of DSMC to include multi-particle collision, in order to make the 
algorithm more efficient in its application to fluids [13]. It is also called stochastic rotational 
dynamics (SRD) [14]. The MPC algorithm consists of alternating streaming and collision 
steps. In the streaming step, the particles move ballistically, r^t + At) = r^(t) + VjAi, where 
At is the time interval between collisions. In the collision step, the particles are sorted into 
cubic cells of lattice constant l c . The collision step consists of a stochastic rotation of the 
relative velocities of each particle in a cell, 

vf w = v c G + H{ Vi - v c G }, (7) 

where v G is the velocity of the center of mass of all particles in the cell. The matrix fl 
rotates velocities by the angle </> around an axis, which is chosen randomly for each cell. The 
translational-momentum and kinetic energy are conserved in the cell. The collision cells are 
randomly shifted before each collision step to ensure Galilean invariance [14]. Although MPC 
fluid originally corresponds to micro-canonical ensemble, the temperature can be controlled 
by an additional rescaling of the relative velocities — v G . In MPC, the angular momentum 
is not conserved and the rotational symmetry is broken by the use of cells. In DSMC [1], two 
particles collide in the cell instead. Thus, the difference between DSMC and MPC is whether 
collisions affect two or all particles in the cell. 

An Andersen-thermostat version of MPC (MPC- AT) has also been proposed [15]. In 
MPC- AT, the velocities in the collision step are obtained as 

vr w = v c G +vr- E v f n /^, (8) 

where vj an is a velocity chosen from a Maxwell-Boltzmann distribution and N c is the number 
of particles in a cell. The summation in Eq. JHJ and the other equations in MPC runs over 
all particles in a cell. The velocity of the center of mass of each cell is conserved, and the 
temperature is constant in MPC- AT (instead of the energy in MPC). 
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Table I - Comparison between DPD and MPC methods of the Langevin dynamics and Andersen- 
thermostat versions. 



DPD MPC 



interacting particles: 






number 


2 


multiple 


chosen by 


distance rij 


collision cell 


potential interaction 


available 


available 


momentum conservation: 






translation 


yes 


yes 


angle 


on/off 


on/off 


energy conservation 


available 


available 



In order to compare the various methods summarized in Table [fl we have to distinguish 
the differences between methods from the variations within each method. Since both of DPD 
and MPC have Andersen-thermostat versions, we compare them first. Several differences 
originate from the variations of each method. The MPC fluid is originally an ideal gas from a 
thermodynamic point of view. It has been generalized by Ihle et al. [17] by a modification of 
the collision rule to produce a nonideal-gas equation of state. However, the equation of state 
can also be changed by the usual potential interactions, when the particles move according 
to Newton's equation mdvi/dt = — Vjf in the streaming step. The angular momentum 
is not conserved in original MPC-AT. However, it can be conserved by the addition of an 
angular-momentum constraint, such that 



V G , v ran _ y ran 



^ ™ n /N c + 1 mil- 1 r J x (Vj - vf n ) \ x r u (9) 

j£cell I j'Ecell 



where II is the moment-of-inertia tensor of the particles in the cell. The energy can also be 



conserved in MPC-AT by the velocity scaling, uf cw -> Uj 2 / XX"$ low ) 2 < cw , where the 

relative velocity is = Vj — vf (MPC-AT-a) or = Vj — vf — {mil J^r, x v 3 } x 
(MPC-AT+a) after the procedure of Eq. JSJ or ©. Note that an additional variable is not 
necessary in the energy-conserved MPC-AT, since undetermined degrees of freedom remain 
for iV c > 3. Thus, two key features can be identified as the genuine differences between DPD 
and MPC: (i) The thermostats act on the relative velocity of two (DPD) or multiple (MPC) 
particles, and (ii) the interacting particles are chosen by their distance (DPD) or by sharing 
the same collision cell (MPC). 

In order to gain a deeper understanding of the relations between DPD and MPC, we 
propose two new methods: DPD with a multibody thermostat (DPD-MT) and MPC-Langevin 
dynamics (MPC-LD). First, we modify the DPD thermostat /dti into a multibody thermostat 
to make DPD more similar to MPC. The number of the thermostats in /dti is 3_/ViV n b/2, 
which is more than the number of degrees of freedom 3N for N n ^ > 2. Since the excess 
thermostats do not play a role, we consider the reduction to three thermostats per particle. 
We define the thermal force in DPD-MT as 
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where w,° = J2j^a w ( r ij)' an d V P — J2j^i w ( r ij) v j / w i 1S the weighted mean velocity. The 
first term of /dt2 is the friction term of /dti, and N n ^/2 thermostats are unified into one 
thermostat between i-th particle and its neighbors. The third and fourth terms are needed to 
conserve the translational momentum. The Fokker-Planck equation for DPD-MT is found to 
be 



dP(X,t) 
dt 



E 



m 



■ di + diTi \ P(X, t) 



T i =w° i {w i -2vl 



knT, 



E 



2d, 



V- w{r jk ) 



dk 



(11) 



where X = {(rj,Vi)|i = l,..,N} and di = d/dvi. The steady state dP(X.,t) /dt = is ob- 
tained in thermal equilibrium. The angular momentum is conserved (DPD-MT+a), when the 
thermostat for the z-th particle is applied only in the direction — rp, where the weighted cen- 



ter of mass is rp = J2j^i y 
can be applied to DPD-MT. 



)rj/w®. We checked that Shardlow's SI splitting algorithm [10] 



Next, we modify the MPC method to the Langevin version (MPC-LD), with 



/MPLT 



-7K-v c G ) + V7Ui(t) 



0(0 



E n ( 

j £ cell 



(12) 



The thermostat is applied to the relative velocities in a collision cell. MPC-LD+a is given by 
the addition of the angular-momentum constraint in the cell, 



/amc 



mU ^ E r J X {iVj-y/liM} 



j/Gccll 



(13) 



The numbers of the thermostats in MPC-LD-a (/mplt) and in MPC-LD+a (/mplt + /amc) 
are 3(7V — iV cc n) and 3(7V — 2iV cc n), respectively, where iV co n is the number of the cells occupied 
by particles. The discrete equation for MPC-LD-a is given by the leapfrog algorithm, 



i(Wi) = vp + a{v 4 (i„) - vp} + b{& 



E 

ji Eccll 



h 



(14) 



where a and b are given by Eq. Eq. (|14fl with r yAt/2m = 1 corresponds to Eq. JSJl of 
MPC- AT. Energy conservation can be added by a rescaling of the relative velocities in each 
cell. Eq. 1|14|) resembles Eq. (0) of MPC. The correlation (vi(t n +i)vi(t n )) decreases with 
increasing r yAt/2m (MPC-LD) or angle <fi (MPC). We checked that the correct (flat) radial 
distribution function of the ideal gas is obtained by all MPC and DPD methods with the SI 
splitting algorithm [10], unlike for some other DPD integrators such as the modified vclocity- 
Verlet algorithm in Ref. [5]. When only two particles are in the cell, the MPC-LD thermostat 
corresponds the pairwise thermostat in DPD. The Langevin versions of DPD and MPC have 
the same relation as the Anderson-thermostat versions. 

The difference between DPD-MT and MPC-LD is the way in which neighboring particles 
are selected. DPD-MT reduces the number of DPD thermostats. However, it does not reduce 
the numerical costs of simulations, since the vectors r^- of all neighboring particles have to 
be calculated. In MPC-LD-a, the numerical costs are reduced; the only information about 
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Fig. 2 - (Color online) Viscosity dependence of MPC-LD-a on (a) the time step At at n = 3 and 
(b), (c) the mean number n of particles per cell at 7* = 1. Symbols represent simulation data, lines 
indicate the results of Eqs. 1)15^ and l|16|l . Viscosity and time are given in units of T70 = y/mk-aT /l c 2 
and to = l c \J mjk-e,T , respectively. The reduced friction constant is 7* = 7T0/TO. Error bars are 
estimated from three independent runs and are much smaller than the size of the symbols. 



particle positions needed is their partitioning into cells. When a very small time step is chosen, 
At <C l c y/m/kBT, particles only move a small distance compared to the cell size in At; in this 
case, the random-shift procedure of collision cells at each time step implies that the particles 
can only react to the time average of the MPC-LD thermostats. This average gives an effective 
weight w m (rij) = \(l — Xij/l c )(l-yij/l c )(l — Zij/l c )\ for \a i3 \ < l c and w sq (r y ) = otherwise, 
where a € {x, y, z}. This is similar to DPD-MT with the weight function w sq (iYj), where the 
rotational symmetry is broken because of w sq (ry ). Thus, MPC-LD can be interpreted as a 
version of DPD-MT, which approximates the weight I0gq(r^) by an uniform weight inside the 
randomly shifted cell. 

Although we have introduced DPD-MT and MPC-LD mainly to fill the missing links 
in Fig. they can be used for practical applications. MPC-LD and MPC-AT need less 
computational costs than DPD for high densities, and have stronger thermostats than MPC. 
As an example, we investigate here the viscosity of MPC-LD-a with the ideal-gas equation 
of state. The viscosities of other methods will be reported elsewhere. The viscosity consists 
of two contribution; the kinetic viscosity 77km and the collision viscosity 77 co i result from the 
momentum transfer due to particle displacements and collisions, respectively. The theoretical 
derivation of the viscosity for MPC [16] can be straightforwardly applied to MPC-LD. The 
viscosities are then obtained as 



7?col 



1 r 



n(l+jAt/2m) 2 At 
(2j/m)(n-l + e- n ) 2 



7(71 - 1 + er n ) 



121 



d-2 



(1 + -fAt/2m) 



(15) 
(16) 



where n = (N c ) and d is the spatial dimension. We also calculate the viscosity from simulations 
for simple shear flow with Lees-Edwards boundary conditions [16, 20] in three dimensions. 
Fig- shows that the theoretical and numerical results are in very good agreement. As At 
or 7/771 increases, 77 kin increases and T7 co i decreases. Thus, the viscosity of MPC-LD can be 
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varied easily. 

In summary, we have proposed new mesoscale simulation techniques — DPD-MT, MPC- 
LD, and variations of MPC-AT — and clarified the relations between several particle-based 
hydrodynamic methods. An obvious question is now which of these methods should be used 
for a given application. The answer depends on the system under investigation and the compu- 
tational demands. For example, the Reynolds number and other dimensionless hydrodynamic 
quantities typically have to be adjusted to match experimental conditions. In general, the 
methods of the MPC group reduce computational costs compared to the DPD group, but 
have the disadvantage of (weakly) breaking the rotational symmetry. We have demonstrated 
here that a comparison of different simulation techniques can stimulate the development of 
new methods, and that ideas developed for one technique can be employed fruitfully in other 
techniques. 

* * * 
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